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ABSTRACT 

Waves observed in the photosphere and chromosphere of sunspots show complex dynamics and 
spatial patterns. The interpretation of high-resolution sunspot wave observations requires modeling of 
three-dimensional non-linear wave propagation and mode transformation in the sunspot upper layers 
in realistic spot model atmospheres. Here we present the first results of such modeling. We have 
developed a 3D non-linear numerical code specially designed to calculate the response of magnetic 
structures in equilibrium to an arbitrary perturbation. The code solves the 3D nonlinear MHD 
equations for perturbations; it is stabilized by hyper-diffusivity terms and is fully parallelized. The 
robustness of the code is demonstrated by a number of standard tests. We analyze several simulations 
of a sunspot perturbed by pulses of different periods at subphotospheric level, from short periods, 
introduced for academic purposes, to longer and realistic periods of three and five minutes. We 
present a detailed description of the three-dimensional mode transformation in a non-trivial sunspot- 
like magnetic field configuration, including the conversion between fast and slow magneto-acoustic 
waves and the Alfven wave, by calculation of the wave energy fluxes. Our main findings are the 
following: (1) the conversion from acoustic to the Alfven mode is only observed if the the driving pulse 
is located out of the sunspot axis, but this conversion is energetically inefficient; (2) as a consequence 
of the cut-off effects and refraction of the fast magneto-acoustic mode, the energy of the evanescent 
waves with periods around 5 minutes remains almost completely below the level /3 = 1; (3) waves with 
frequencies above the cut-off propagate field-aligned to the chromosphere and their power becomes 
dominating over that of evanescent 5-minute oscillations, in agreement with observations. 
Subject headings: MHD; Sun: chromosphere; Sun: oscillations; Sun: photosphere; sunspots 



1. INTRODUCTION 

Sunspots can be considered as laboratories for studies 
of magnetized plasma in conditions that are inaccessi- 

ON ble on Earth. They give clues about the physics of en- 
ergy propagation (e.g., in the form of different oscillatory 
modes) in fluids permeated by strong magnetic fields. 

\q o From the analysis of wave propagation it is possible to 
infer basic properties of stellar atmospheres. Waves ob- 
served in active regions are believed to have a relevant 
. role in the energy balance of the solar atmosphere, being 
one of the candidates to explain the chromospheric heat- 

. £^ ' ing . The first detecti o n of w aves in sunspots was done 
bv iBeckers fc Tallantl (|1969[ ). who called them "umbral 
flashes" . Since this date waves have been found in many 
solar features, from photospheric flux tubes to the solar 
wind. 

Oscillations show a different behavior in different re- 
gions of sunspots. In the umbral photosphere their 
power spectrum peaks around 3 mHz (period of 5 min- 
utes), and their amplitudes are around a hundred m s _1 . 
Waves in the photosphere of sunspot umbrae are sim- 
ilar to those observed in the qui et sun, but with thei r 
power reduced by a factor of 2-5 (Woo dsfc Craml fl981). 
Several mechanisms have been proposed to explain this 
power suppression: reduction of wave excitation inside 
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sunsp ots (Go ldreich fc Keeleylll977tlGoldreich fc Kumar 
1988, ll99Qh . ff-mode absorption inside sunspots (|Callv 
1995), different heights of spectral line formation due to 
the Wilson depression and altering of p-mod e eigenfunc- 
tions by the magnetic field (Hi ndman et al.lll997[ ). 

Chromospheric umbral oscillations have a major power 
peak around 5 mHz (period of 3 minutes) and amplitudes 
of several kilometers per second. Velocities measured 
in chromospheric umbrae sho w saw-tooth temporal pro- 
files, typical for shock waves (Lites 1986; Centeno et al.l 
2006). To explain the multiple peaks in the 3 minute 
band in the power spectrum of the chromospheric waves 
several hypotheses ha ve been investigated: a resonant 
chromospheric cavity (Zhugzh da et al .1119851): non- linear 
inter action of photospheric modes (jGurman fc Leibacherl 
1984); slow magneto-acoustic mode field-aligned propa- 
gation from the photosphere to the chromosphere in the 
5-6 mHz band (jCenteno et al.ll2006D . 

In the penumbra the dominant oscillatory phenomena 
are running penumbral waves, which are best seen in Ha 
core velocity observations as disturbances propagating 
radially outward from the umbra. As they move across 
the penumbra their radial velocity is apparently reduced 
and their frequencies also decrease from 4 — 5 mHz near 
the umbral/penumbral boundary to aro und 0.7—1.5 mHz 
at the outer edge of the penumbra (Bogdan fc Judge! 
l2006h . 

It is now becoming apparent that the photospheric 5- 
minute oscillations, chromospheric 3-minute oscillations, 
and running penumbral waves are the different man- 
ifestation of the same global dynamical phenomenon 
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(jvan der Voort et al.ll200l iBloomfield et all 12007( 1. No 
fully comprehensive model of this global oscillation phe- 
nomenon was presented up to date due to the compli- 
cated mathematical description of the physical processes 
playing a role in realistic magneto- atmospheres. The re- 
cent progress in observations and numerical simula tions 
of sunspot waves is summarized in iKhomen ko1 (|20T0h . 

Although a lot of a nalytical work was done 
in simple atmospheres (|Ferraro fc Plumptonl |l958; 
Zhugzhda & Dzhalilov 1984, to name a few), those 
works were restricted to very idealized cases. Nu- 
merical simulations allow more flexibility and over the 
last years many attempts were done to perform nu- 
merical modelling of waves in non-trivial magnetic 
field configurations, with applicati ons to the photo- 
sphere and the chromosphere (e.g . [Car gill et al.l 119971: 
Rosenthal et al.l 120021: lHasan e t al. 2003; Bogdan et al.l 
20031: lHasan fc U lmschneider 2004; Has an et al.l 120051 : 
Khomenko fc Colladosl 12000 : iKhomenko et al]l2008jT ln 
most of the cases, these works were restricted to stud- 
ies of the behavior of high-frequency waves (above cut- 
off) in two-dimensional situations. Despite these limi- 
tations, several important questions were learned from 
these models. It was shown that the fast magneto- 
acoustic mode in the magnetically dominated region (i.e. 
where the sound speed cs is much lower than the Alfven 
speed va) is refracted and reflected down to the gas 
pressure dominat ed atmosphere due to the gr adients of 
the Alfven sp eed (Khome nko fc Collad os 2006). Earlier, 
iRosenthal et al.l (j2002f ) demonstrated that the inclina- 
tion of the magnetic field lines is important for the fast 
mode reflection, i.e. in those regions where the incli- 
nation angle is large almost all of the fast mode wave 
energy is reflected back down. Another important fea- 
ture of all these simulations is the presence of the mode 
transformation at the layer where va = cs- Around 
this layer, the phase speeds o f all modes are similar and 
different waves can interact dBogdan et al.l l200l ICalivl 
2006; IKhomenko fc Collados! |2006|). The direction and 
efficiency of the mode transformation depends on the fre- 
quency of the wave and th e angle between t he wave vec- 
tor and the magnetic field (Cally 2005, 2006). When this 
angle is small and the frequency is high, the fast mode 
can be converted into the slow mode and vice versa. 

Most of the works cited above on the simulations 
of wave propagation in the upper layers of sunspot s, 
as well as the analytical theories of the mode trans- 
formation, were developed for high-frequency waves 
with frequencies above, or just at, the cut-off fre- 
quency. There is another class of simulations where 
the problem of helioseismic wave propagation below 
sunspots is addressed and waves with realistic so- 
lar fr equencies in the 3-5 mHz range are consid- 
ered dCallv fc Bogdanl 119971: IParchevskv fc Kosovichev 
l2009HHanaso ge 2008; Camero net al.ll2008l : iMoradi et al. 
2009 rj IKhomenko et al.l l2009). In the latter works most 
attention has been paid to the wave propagation in sub- 
surface layers, rather than in upper photospheric and 
chromospheric layers. Thus, there is gap between these 
two kind of models. With the development of our code 
and its first results presented here we pretend to cover 
this gap and to address the problem of three dimensional 
propagation and transformation of the 3-5 mHz period 
waves in the upper layers of sunspots. Our aim is to 



investigate the response of the magnetic atmospheres to 
oscillations with low frequencies (like 5-minute waves) 
and to compare this response with that from the high- 
frequency perturbations. As those 5-minute waves have 
most power at the photospheric heights and are used in 
measurements to derive the parameters of the solar at- 
mosphere, including solar active regions, understanding 
their propagation and transformation properties is im- 
portant. Our objective is to perform numerical calcula- 
tions sufficiently realistic to imitate the wave excitation 
in sunspots, reproduce the change of wave frequency with 
height, formation of shocks at chromospheric layers, and 
to be at the level allowing comparison with photospheric 
and chromospheric observations by spectral synthesis. 

Several numerical codes have been developed among 
the solar physics community to study the propagation 
of waves in magnetized atmospheres. They adopt differ- 
ent strategies for the numerical scheme, boundary con- 
ditions and w ave driving, having t heir advantages and 
disadvantages (jMoradi et al.l 1200 9 at ). The upper bound- 
ary typically represents a problem in wave simulations, 
since waves should not artifi cially be reflected there back 
into the physical domain. IRosenthal et al.l ([2002) ap- 
ply characteristic boun dary conditions at the top bound- 
ary; [Hasan et al. (2005) uses t he open boundary concept; 
while IKhomenko et al.l ([2008) introduced a special me- 
dia at the top called Perfectly Matched Layer (PML) 
which absorbs with almost no reflections the waves that 
reach the upper bo undary. IRosenthal et al.l ([2002) and 
lHasan et al.l (|2005|) solve the complete MHD equations, 
while IKhomenko et al.l (j2008h solve equations for pertur- 
bations with all non-linear terms retained. This strategy 
gives them an advantage for precision of the numerical 
scheme and for the application of the boundary condi- 
tions. 

3D MHD codes for wave s imula tions are starting to 
be available. Ca meron et al.l (|2007| ) presented the semi- 
spectral linear MHD code SLiM, developed for helioseis- 
mology purposes. In this code, the horizontal derivatives 
are evaluated in Fourier space and the upper bound- 
ary is treated as a sponge layer. Another 3D linear 
MHD code for wave propagation has been developed by 
IParchevskv fc Kosovichevl (|2QQ7h . The authors use the 
realistic OPAL equation of state and a PML layer as the 
upper boundary condition. Numerically, the upper mag- 
netized atmospheric layers represent an additional prob- 
lem limiting strongly the time step of the simulations 
due to the high values of the Alfven speed. To overcome 
this problem, one of the strategies used is the Lorentz 
force controller. This method consists in reducing the 
amplitude of the Lorentz force in the layers where the 
Alfven speed is large. This method is used by iHanasogd 
(2008) in his 3D linear code. However, the influence of 
this artificial procedure on the simu lated wave properties 
has not been verified yet. Recently Shelva g" et al.l ([2008) 
presented a nonlinear 3D parallel code developed on the 
base of the VAC code (Toth 1996). In this code, they use 
the same philosophy as IKhomenko et all {2008), solving 
non-linear equations for perturbations. As for sunspots, 
only results of the 2D helioseismic wave propagation be- 
l ow surface have be en presented so far from this code 
(Shel vag et al.ll2009l h not addressing the wave propaga- 
tion in the upper layers. 

The non-linear MHD code described in this paper rep- 
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resen ts an extension of the code by iKhomenko et al.l 
(2008) into three dimensions. Our code solves the MHD 
equations for perturbations in a background plasma in 
magnetohydrostatic equilibrium using artificial diffusion 
to assure numerical stability of the solution, and is fully 
MPI-parallelized. In this paper we present extensive nu- 
merical tests demonstrating the code stability and ro- 
bustness. We also present the first scientific results of 
the three-dimensional wave propagation and transforma- 
tions in the upper atmosphere of a sunspot model. Our 
first aims are to study the propagation of waves, excited 
by sub-photospheric pulses, up to the chromosphere, in- 
cluding the shock wave formation. We present a detailed 
description of the three-dimensional mode transforma- 
tion for waves with different frequencies, from below to 
above the acoustic cut-off frequency. In particular, we 
address the question of conversion to Alfven waves in re- 
alistic conditions. Works on conversion from/ to Alfven 
wav es in three dimensions a re only now being initiated 
(e.g. iCallv fc Goos sens 2008) and no generalized picture 
has been obtained yet. The organization of the paper is 
as follows. In Section [2] we describe the MHD equations 
and the numerical details of the calculations. In Section 
02 we discuss the results of the simulations of the three- 
dimensional wave propagation and mode transformation 
in a sunspot model. This section includes the descrip- 
tion of the magneto-static sunspot model in equilibrium 
( Section f3TTj) and a brief discussion of the di sper sion rela- 
tions and wave propagation speeds (Section 1X2]) . Section 
[H presents our conclusions and our future plans for the 
application of the code. Finally, in Appendix [Aj we show 
several numerical tests applied to prove the robustness of 
the code. 

2. NUMERICAL PROCEDURE 

2.1. MHD equations 

Our numerical code is an extension to 3 dimension s 
of the code desc r ibed i n IKhomenko fc Colladosl ((2006); 
IKhomenko et al.l (j2008[ ). It solves the non-linear equa- 
tions of ideal compressible MHD. Written in conservative 
form, these equations are: 
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dt 
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P(g-V) + Qrad, 

(3) 
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where I is the identity tensor, p is the density, v is the 
velocity, p is the gas pressure, B is the magnetic field, g 
is the gravitational acceleration and e is the total energy 
per unit volume, 
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(5) 



The dot '•' represents the scalar product of vectors, 
while the notation 'BB' stands for the tensor product. 
The energy losses Q r ad can be described by the Newton 



law of cooling, but, in the simulations presented in this 
paper, it is set to zero. We neglect the viscous force, the 
thermal conduction and the terms describing the diffu- 
sion of the magnetic field. However, artificial equivalents 
of some of these terms are introduced later for the issue 
of numerical stability of the simulations. The term S(t) 
in Eq. [2] represents a time-dependent external force. 

In an equilibrium state, where temporal derivatives 
and velocities are null, and in the absence of external 
forces (S = 0), the previous equations reduce to the equa- 
tions of the force balance for a gravitationally stratified 
magnetized plasma: 
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[Po 



2/W 



BqBq 
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(6) 



Considering departures from the equilibrium state in- 
duced by an external force S, p, p and B can be expressed 
as the sum of the background value (subindex 0) and the 
perturbation (subindex 1): 



and 



Po + pi , 



P = Po + Pi , 



B = Bq + Bi , 



(7) 



(8) 



(9) 



while the velocity only corresponds to a perturbed value 
v = Vi. 

The non-linear equations for perturbations are ob- 
tained by replacing expressions ([7H9]) in Eqs. (p~H!]) and 
subtracting the equation of the magnetohydrostatic equi- 
librium (Eq. [6]). The following system of MHD equations 
for perturbations of density, pressure, magnetic field and 
velocities is obtained in conservative form: 
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dt VX [ VlX ( Bo+B ^ + (f 1 J dlff ' < 13 ) 

Artificial diffusion terms have been added to Eqs. (fTTI - 
[T3|) compared to Eqs. (j2]-|4]). The diffusivity terms in 
Eqs. (pTHT3|) have their physical counterparts and are 
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needed for reasons of stability of the simulat ions. A simi- 
lar st rategy is applied in the MURAM code (jVogler et al.l 



The code solves the above system of non-linear equa- 
tions for perturbations. The use of equations for per- 
turbations instead of the complete equations for wave 
simulations has two big advantages. Firstly, the terms 
describing the static model and those for perturbations 
can vary by orders of magnitude. Thus, by excluding 
equilibrium terms we avoid important numerical preci- 
sion problems. Secondly, the boundary conditions are 
e asier to implement on equations for perturbations (see 
§& 

In high layers, where the magnetic pressure is much 
larger than the gas pressure, Eq. ([T2]) is numerically 
problematic as recovering thermal energy (p) from total 
energy (e) leads to numerical errors. In these layers it 
is replaced by the equation describing the balance of the 
internal energy 



dpi 
dt 



• viVOo +4 v((>o + Pi)vi 

( dpi 



-ViV(po+Pl) =(7-l)PlQrad+( : ^) . .(14) 

The transition from one equation to the other is done 
smoothly using the plasma parameter f3 as a criterion. 

The gravity g is constant with height and the equation 
of state of an ideal gas with constant 7 is used to link 
the thermal variables. 

The computational domain is discretized using a three- 
dimensional Cartesian grid with constant space step in 
each dimension. The spatial derivatives are approxi- 
mated by a centered, fourth order accurate, explicit finite 
differ ences scheme using five grid points (|Vogler et al.l 
2005). The solution is advanced in time by an explicit 
fourth-order Runge-Kutta. The fourth order differences 
allow to give an accurate solution, still resolving shock 
fronts, which is difficult with higher-order methods. 



The time step has to ensure that the physical depen- 
dence domain is included inside the numerical depen- 
dence domain. According to this, the mesh width must 
be bigger than the distance traveled by the information 
in a single time step due to mass flow, waves or diffusion 
transport. The time step must be chosen to be smaller 
than the advective time step and the time step imposed 
by the diffusion terms: 



At < min(At v , At diff ) 



(16) 



In this expression At v is the time step imposed by a 
modified CFL condition, approximately valid for MHD 
equations, 



At v 



1/Ax 2 + 1/Ay 2 + 1/Az 2 



1/2 1 



(17) 



where v max is the maximum value of the sound and 
Alfven speeds. The time step imposed by the diffusion 
Atdiff corresponds to the minimum of the diffusion time 
across the three dimensions, 



Atdiff = Cdiff min ( 



Ax 2 Ay 2 Az 2 



(18) 



where the constant coefficients c v and Cdiff are taken to 
be below one to ensure the stability of the solution. The 



diffusion coefficients v x ,y,z are those defined in § 12.21 
2.4. Filtering 

In the particular case of wave simulations high diffusion 
is not desirable since it modifies the wave amplitudes. At 
the same time, low diffusion can not always prevent the 
developing of high frequency noise. In such cases we per- 
for m an additional filtering of s mall w avelengths. Follow- 
ing (Parchevsk^^Sosovich^ (|2007f ) we use a sixth-order 
digital filter to eliminate unresolved short-wave compo- 
nents: 



2.2. Artificial diffusivity 

To damp high-frequency numerical noise on sub grid 
scales, we replace the physical diffusive terms in the equa- 
tions of momentum and energy by artificial equivalents. 
In the induction equation we retain the magnetic dif- 
fusion term, replacing n by an artific ial value. In gen- 
eral, we use a philosophy sim ilar to iStein fc Nordlundl 
(fl998h . ICaunt fc Kor^il (|200lh and iVogler et al.l (|2005h . 
Each physical quantity has its own diffusivity coefficient 
(scalar/vectorial for scalar /vectorial quantities), which is 
formed by a shock resolving term, a hyperdiffusivity part 
and a constant contribution: 



n {u) = vt\u) + v^{u) + vl 



(15) 



where uf = (cs + va)AxiF(x, y,z), u is the correspond- 
ing quantity, F(x,y,z) gives the form of the constant 
contribution and Axi = Ax, Ay, Az. 

2.3. Time step 



umt = u(x) 



E 



d m u(x + mAx), 



(19) 



where u is a variable before filtering and umt is after 
filtering. The filter can be applied in the three spatial 
coordinates independently. The coefficients d m have been 
chosen to construct the filtering function: 



G(kAx) = 1 - J2 dme imkAx = 1 - sin ( 



m— — 3 



/kAx 



. (20) 



The frequency of application of the filter depends on the 
simulation, but it is usually applied every 10 seconds. 

2.5. Boundary conditions 

Boundary conditions are an important issue for wave 
simulations. One usually wants to prevent spurious wave 
reflections at the boundaries. Two strategies commonly 
applied are based on characteristic boundary conditions 
or sponge layer. Ca lculating characteristic conditions 
(Ros enthal et al.ll2002l ), apart from tricky, gives good re- 
sults in simple magnetic field configurations, when the 
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wave propagation directions are easily predictable. For 
more complex magnetic field configurations, the calcula- 
tion of the characteristic directions is not an easy task. 
The other alternative, i.e., sponge layer, consists in lo- 
cating an absorbing layer at the boundary to dissipate 
the wave energy and prevent it from coming back to the 
physical dom ain. This strategy is implemented in the 
code SLiM (jCameron et al.ll2007[ ). Absorbing layers give 
goods results only when the absorption is gradual and 
needs a large amount of grid points. Thus, numerically 
they are very costly. In our code we use d yet another a l- 
ternative, the Perfectly Matched Layer (jBerengerl LL994). 

The Perfect Matched Layer (PML) is designed to ab- 
sorb waves without reflection s. This method was first in- 
troduced by [Berengerj (|1994l ) to absorb electromagnetic 
waves in numerical solutions of Maxwell equation s. Later 
it has been applied to Euler equations (Hu 1996) and to 
simulations of acousti c waves in a strongly stratified so- 
lar convection zone (jParchevskv fc Kosovic hev 2007|). In 
our code we extend the method to the full set of the 
MHD equations (Eqs. fT0HT3j) . The MHD equations can 
be written schematically, in conservative form, as: 



du <9F(u) <9G(u) <9K(u) 



dt dx 



dy 



dz 



H(u), (21) 



where u = [pi, (po + pi)vi, ei, Bi] is the vector that con- 
tains the conserved variables; the vectors F, G and K are 
the fluxes, whose expressions can be found in the system 
([TQUT3]) : and H represents the source terms at right-hand 
side of the same equations. Inside the PML, variables u 
are split into three components in such a way that u = 
ui + u 2 + u 3 and also H(u) = Hi(u) + H 2 (u) + H 3 (u) 
and the system of MHD equations is split into a set of 
three coupled, one dimensional equations: 



■<7 x (a0ui =Hi(u), (22) 



<9ui 


l_ dF(u) 




dx 


du 2 


dG(u) 




dy 


<9u 3 


3K(u) 


dt 


dz 



■<Ty(y)u2 = H 2 (u), (23) 



a^)u 3 = H 3 (u). (24) 



These split Eqs. (|22H2^|) are solved independently in 
the PML, in contrast to the unsplit forms (Eq. |2T|) which 
are solved in the physical domain. 

An absorption term has been added to each equation. 
The coefficients <j x (x) /<j y (y) /a z (z) only depend on x/y/z 
coordinate and are non-zero in the x/y/z PML faces, re- 
spectively. Theoretically, a PML with constant absorp- 
tion coefficient produces no reflections for plane waves 
incident on a flat interface for any angle of incidence 
and any frequency. However, due to the finite differ- 
ence implementation of the PML equations in numerical 
calculatio ns, reflect i ons m ay appear when a has a steep 
gradient (Be rengerl LL996). To solve this problem it is 
necessary to include smooth variations in the absorption 
coefficients from small values at the interface between the 
PML medium and the physical domain to large values at 
the outer boundary. Following iHul (|2001[ ), good results 
are obtained with absorption coefficients of the form: 



a / 
Ax V 



X - XpML V 
^PML ' 



b r y - 2/pml V 
Ay \ ypuL J 

c / z-zpuL y 

Az V Zpy[L ' 



(25) 



(26) 



(27) 



where Ax, Ay and Az are the discretization steps, a, b 
and c are constants controlling the damping amplitude, 
and xpml, i/pml and zpml are the thickness of the PML 
domain in each spatial direction. In a typical calcula- 
tion we need a PML with 10-15 grid points. The coef- 
ficients a, b and c depend on each particular simulation 
and are proportional to the wave speed at the boundaries. 
We locate PMLs at all boundar ies of our simulat i on do - 
main. The results presented in iKhomenko et al.l ([2008) 
show that this strategy gives good results even for strong 
shocks. 

The PML formulation, as presented above, can be- 
come unstable in long simulation run s. According to the 
literature (see, e.g., Hest havenlll998l ). a high frequency 
noise filtering can improve the long-time stability of the 
PML layer. We found that this method works well for 
MHD waves. Applying a filter allows to delay the ef- 
fects of the possible instability the necessary time to 
complete long enough simulation runs. Previously the 
PML layer for MHD wave simu lations was applied by 
iParchevskv fc Kosovichevl (|2009f h though the stability of 
the PML formulation for MHD equations was not dis- 
cussed in this paper. 

The code can also account for periodic and closed 
boundary conditions. 

2.6. Parallelization 

Parallelization has been done with MPI following a dis- 
tributed memory concept in which all data used by a 
processor are situated on the memory partition accessi- 
ble to it. Data are split in a certain number of processors 
by means of a domain decomposition scheme. The full 
numerical domain is divided into a set of three dimen- 
sional subdomains, with communication between proces- 
sors only occurring at their common data boundaries. 
For this purpose, each domain includes three layers of 
"ghost" cells at each boundary. The 5-point stencil of the 
fourth-order scheme needs two cells outside the subdo- 
main, while for the filtering it is necessary to include one 
more "ghost" layer. The boundaries of the subdomains 
which are neighbors of other subdomains receive directly 
the required information and store it in the "ghost" lay- 
ers, while when the boundary of the subdomain coincides 
with the global boundary, the "ghost" layers are settled 
with the values imposed by the boundary condition. 

We have performed a number of standard numerical 
tests, as well as tests specifically devised for wave propa- 
gation, to prove the robustness of the numerical method 
and of the boundary conditions. The results of these 
tests are given in Appendix [A[ 



3. 



3D PROPAGATION AND TRANSFORMATION OF 
WAVES IN SUNSPOT MODEL 



Below in this Section we discuss several simulations of 
the propagation and transformation of MHD waves in a 
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magnetostatic sunspot model, excited by pulses with dif- 
ferent periods and locations. To facilitate the reading of 
this section, Table 1 summarizes the simulation runs. It 
gives the number of the sub-section where the results are 
presented, the properties of the driver, its location, hor- 
izontal source size, the duration of the simulations and 
the numbers of the corresponding figures. In all cases, 
to identify the different wave modes in three spatial di- 
mensions we use projections of the velocity into three 
characteristic directions. To quantify the mode transfor- 
mation we calcula te t he acoustic and magnetic energy 
fluxes (see Section [373]) . 

We use a vertical force in the momentum equation S(t) 
to perturb a magnetostatic sunspot atmosphere in equi- 
librium and study the waves generated by this perturba- 
tion. We have performed several numerical simulations, 
all of them with the source situated below the quiet pho- 
tosphere at z = —0.5 Mm and y = Mm, but with 
differences in the period, the horizontal x location of the 
source relative to the axis of the sunspot, as well as the 
horizontal size of the source R svc . 

In the simulations described in Sect. 13. 4[ 13.51 and 13.61 
the temporal behavior of the driver is harmonic is de- 
scribed by the expression: 

lirt 

S z (r, t) = AP(r) sin (28) 

r 

In this equation, A is the amplitude of the 

source, P(r) = [l — (r/R src ) 2 ] describes the 
source shape, R SYC is the source radius, r = 
y/(x- x src ) 2 + (y - Vsvc) 2 + (z - z src ) 2 is the dis- 
tance from the source center and r is the period of the 
harmonic source. P(r) is zero if r > R S rc- The x and y 
components of S, S x (r,t) and S y (r,t), are set to zero. 

In the simulations described in Section [3771 the behavior 
of the driver is not harmonic in time, but rather has a 
shape of Ricker wavelet: 

S z (r, t) = AP(r) (1 - 2r 2 ) e~ T ° (29) 

where To = uJot/2 — tt. Such driver produces 
a spectrum of waves with a central frequency ujq 
(jParchevskv fc Kosovichevl I2009D . We set uj = 3.33 
mHz, so the spectrum of our driver resembles a solar 
one and covers a broad range of frequencies. 

This latter run is particularly interesting, as it allows 
us to study the behavior of a realistic spectrum of so- 
lar waves in the upper layers (photosphere and chromo- 
sphere) of a sunspot model, including the propagation 
of individual wave modes and wave energy fluxes. As 
far as we are aware of, no such investigation has been 
performed as of today. 

3.1. Magnetostatic sunspot model 

Here we use a magnetostatic (MHS) model atmosphere 
in equilibrium represen t ative o f a sunspot, adopted from 
iKhomenko fc Colladosl ((2008). This MHS model is 
a thick flux tube with distributed currents, it is az- 
imuthally symmetric and has no twist. The variations 
of field strength and gas pressure are continuous across 
the spot. At 40 Mm far from the sunspot axis the model 
merges smoot hly into a quiet Sun atmosphere ta ken from 
the model S (iChristensen-Dalsgaard et allll996[ ) in the 



deep sub- photosphere layers a nd continuing as a VAL- 
C model (jVernazza et al.|[T981h in the photospheric and 
chromospheric layers. The sunspot axis in the atmo- 
spheric layer s is given by the semi-empirical model of 
Avrett (1981). The dimensions of our computational do- 
main are 2.5 Mm in the vertical direction and 15 Mm in 
each horizontal direction with a grid size of Az = 25 km 
and Ax = Ay = 75 km. The bottom level is 1.25 Mm be- 
low T5000 = 1 at the photospheric quiet sun atmosphere, 
which was chosen as the zero level of the coordinate z. 
The top level is 1.25 Mm above z = 0. A PML layer of 10 
grid points was used at both bottom and top boundaries. 
With this, the physical domain occupies from z = — 1 
Mm to z = 1 Mm. The axis of the sunspot is placed at 
the center of the simulation domain. The magnetic field 
at the axis is about 900 G at z — Mm. 

3.2. Dispersion relations and propagation directions 

In a three-dimensional situation like the one consid- 
ered here, three types of wave modes exist: fast and slow 
magneto- acoustic waves and the Alfven wave. Each of 
these modes is described by its own dispersion relation, 
i.e. relation between its temporal frequency and the wave 
number uj = cj(k). The modes are characterized by their 
own phase velocity, v p h = tj/k, and the group velocity, 
v g = <9cj/<9k. The first one is parallel to the direction of 
k and gives the wave front propagation velocity. The sec- 
ond one defines the direction of the energy propagation. 
The magnitudes and directions of v p h and v g usually do 
not coincide. 

In a 3D homogeneous unstratified atmosphere perme- 
ated by a constant magnetic field B , parallel to z axis, 
and y axis normal the plane defined by the direction of 
propagation k and Bo, the generalized wave equation for 
the perturbations decouples into two, defining the disper- 
sion relation for the three MHD modes: 

%■ = Va COS if, (30) 

k 

^ = \(v\ + 4)±\y/(v% +C|)2 -AV\C% COSV (31) 

where va = (Bq/ l^opo) 1 ^ 2 is the Alfven speed, cs = 
(TPo/po) 1 / 2 is the sound speed, and tp is the angle be- 
tween k and Bo- 

Equation ([30]) describes an Alfven wave, whose associ- 
ated perturbations in vi and Bi are transversal to Bq 
and k and the perturbations in pressure and density are 
zero. For the propagation along the magnetic field the 
phase velocity of this wave is equal to va- The energy of 
the Alfven mode propagates along magnetic field lines, 
according to its group velocity. 

The solutions with the plus and minus sign in Eq. ([3T]) 
are the fast and slow modes, respectively. When one of 
the cs or va is much higher than the other, the dispersion 
equation (Eq. [31]) can be simplified to uj = kvf ast for the 
fast mode and uj = kv s \ ow cos ip = v s iowk • Bo/|i?o | f° r the 
slow mode, where Vf ast is the highest velocity between 
va and cs and v s \ ow is the lowest. The direction of the 
phase and group velocities of the fast mode is k and their 
magnitude is either cs of va, so the energy propagates 
in the same direction as the wave front. The slow wave 
group velocity is directed along Bo- When cs and va are 
similar, the behavior of the phase and group velocities of 
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TABLE 1 

Summary of the simulation runs 



Section 


Driving 


Location 


Rsrc 


Duration 


Figures 


4.4 


50 s, harmonic 


x = Mm 


150 km 


904 s 


7-12 


4.5 


50 s, harmonic 


x = -3 Mm 


150 km 


1023 s 


13-17 


4.6 


180 s, harmonic 


x = Mm 


540 km 


930 s 


18 


4.7 


300 s, wavelet 


.x = Mm 


900 km 


1511 s 


19-21 


4.7 


300 s, harmonic 


.x = Mm 


900 km 


1212 s 





the fast and slow modes deviates from this simple picture, 
and can be retrieved from the derivation of v g from Eq. 

SB. 

3.3. Identification of the wave modes in simulations 

In the case of a real atmosphere, the division into pure 
wave modes is not so simple as in the idealized case 
described above, as often no clear distinction between 
the modes can be done neither physically nor mathe- 
matically. Even in the above simple case the governing 
partial differential equation factors into a single second 
order wave equation for the Alfven mode and a fourth 
order wave equation for the coupled fast-slow modes, so 
the idea that there are three distinct modes may not al- 
ways be correct. However, the simplicity of this picture 
makes it attractive and we will discuss the properties of 
the waves in realistic atmospheres in terms of the three 
modes. 

To help the identification of the wave modes in simu- 
lations, we use the mode properties described above. We 
project the vectorial quantities (velocity and magnetic 
field perturbations) into the directions aligned/normal 
to the equilibrium magnetic field B . 

At each location of the computational domain, we cal- 
culated the projections of the vi and Bi into the follow- 
ing Cartesian directions: 

eiong = [cos (j) sin 0, sin sin 0, cos#], (32) 

e P erp = [ — cos (j) sin 2 sin </>, 1 — sin 2 sin 2 </>, (33) 
— cos sin sin 0] , 

etrans = [ — cos 0, 0, cos (J) sin 0] . (34) 

where is the magnetic field inclination from the verti- 
cal and <p is the field azimuth, measured from the x — z 
plane. The direction of ei on g is along the magnetic field 
Bo- The directi on of e ner r> is normal to th e field and was 
chosen following iCallv fc Goos sens (20Q8|) as the asymp- 
totic polarization direction of the Alfven mode in the 
low-/3 regime. The last component etmns = ei ong x e pe rp 
is set in the direction normal to the other two. 

We expect that in a region where cs > va, the slow 
magneto-acoustic mode will be identified in e tra ns pro- 
jection of the velocity vector, while the fast magneto- 
acoustic mode will be equally visible in all velocity com- 
ponents as is propagates isotropically. In a region where 
cs < va-> the slow magneto- acoustic mode will be identi- 
fied projected into ei ong direction, the Alfven mode pro- 
jected into the e pe rp direction, and the fast magneto- 
acoustic mode in the the direction normal to these two. 



To quantify the amount of energy contained in differ- 
ent wave modes and to develop a measure of the mode 
transformation, suitable in the case of complex mag- 
netic field configurations like the one considered here, 
we found it useful to cal culate the wave energy fluxes 
(Bray & Loughhead 1974). The acoustic energy flux is 
given by the expression: 

F ac =PiVi, (35) 

and magnetic energy flux is given by: 

F mag = Bi x (vi x B )//i . (36) 

The acoustic energy flux contains the energy of the 
wave with acoustic nature, which corresponds to the fast 
mode in the region where va < C5, and to the slow mode 
in the region where va > cs. In this region, the magnetic 
flux includes the fast and the Alfven modes. Since, as we 
will see in the next section, in the region above the layer 
where va = cs the fast mode is refracted down towards 
the photosphere, the magnetic energy which propagates 
upwards along field lines must correspond to the Alfven 
wave, making possible the identification of this mode. 

3.4. Case of 50 s harmonic force located at the axis 

Figure [T] presents a two-dimensional snapshot of some 
variables and Figure [2] gives the temporal evolution of 
the projected velocities in the simulation run with the 
harmonic 50 sec force located at the sunspot axis. The 
panels (a-c) in Figs. [I] show the longitudinal, transversal 
and perpendicular velocities scaled with a factor ^/pQ > . 
These magnitudes provide the square root of the kinetic 
energy associated to the waves. Due to the strong density 
fall off, some velocity perturbations at high layers have 
so low energy that makes them indistinguishable in this 
representation. To make them visible, we have plotted 
additional contours of constant velocity. In the case of 
the pressure (Fig. []J), its drop with height makes the 
absolute valued of the perturbations at the lower layers 
much higher than at the upper layers. In this case, the 
contours represent the ratio of constant pi/po. 

3.4.1. Propagation below the surface 

The vertical force acts in a region where c 2 s /v\ ~ 9.1 
and it generates mainly an acoustic fast mode, whose 
oscillations can be seen in the longitudinal velocity, pres- 
sure and density snapshots in Fig. [TJ This vertical im- 
pulse produces initially a deficit in density and pressure 
at the place where the source is located and, because of 
that, horizontal motions also appear, creating a magnetic 
slow mode seen in the transversal velocity and magnetic 
field snapshots in Fig. [TJ 

At photospheric level, the longitudinal velocity has an 
amplitude of about 200 m s _1 , the amplitude of transver- 
sal velocity is 50 m s _1 and the transversal magnetic field 
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RADIAL DISTANCE [Mm] 

Fig. 1. — Variations of the velocity in the direction ei on g (a), etrans (b), and e per p (c), all of them scaled with factor ^/po of the unperturbed 
density; magnetic field in the direction etrans (d), and pressure (e) at an elapsed time t = 820 s after the beginning of the simulations for 
the 50 s harmonic force located at x = km, y = km and z = —500 km. All panels show the plane y = km, except panel (c), which 
shows y = 400 km. Black inclined lines are magnetic field lines. Green lines are contours of constant v^/cg- The image color coding is 
such that blue colors represent lower values and red colors are higher values with respect to the mean. Scaling in panels (a), (b), (c) is the 
same. Figure (a) shows contours of equal longitudinal velocity, and (d) shows contours of equal pi/po- 



oscillates with a maximum deviation from the equilib- 
rium value of 4 G. 

The temporal evolution given in Fig. [2] shows the fast 
mode (acoustic in nature) propagating in the deep lay- 
ers upwards to the region where va ~ cs- It appears as 
a perturbation in the longitudinal velocity (panels c-d). 
Variations in longitudinal velocity are accompanied by 
acoustic variations of pressure and density (not shown in 
the figure). The slow mode (of magnetic nature) is vis- 
ible in the transverse the acoustic perturbation reaches 
the surface v\ = c| earlier than the magnetic perturba- 
tion. In these deep layers, the acoustic oscillations have 
a wavelength larger than the magnetic ones. As the fast 
(acoustic) wave propagates in the region va < cs, its en- 
ergy is distributed in the three spatial dimensions and it 
decreases away from the source as 1/r 2 . Once it reaches 
the va > cs region, the energy redistribution in horizon- 
tal directions is not so important because it is channeled 
along the field lines and is only affected by the density 



falloff. 



3.4.2. Three-dimensional mode transformation 



When the waves reach the v\ 



c| layer from below, 
several mode transformations take place in the simula- 
tion. 

First of all, the fast acoustic mode moves from a re- 
gion where va < cs to another where va > cs keeping 
its acoustic nature but changing from fast to slow mode. 
This transformation can be seen in the snapshots of lon- 
gitudinal velocity and relative pressure in Fig. [T] (panels 
a and e) as the wavefronts above the layer v\ = c|. The 
atmosphere above v\ = c| is dominated by the mag- 
netic field and this slow acoustic mode propagates up- 
wards along field lines (Fig. [2j a-c). The amplitude of 
this wave increases according to the density drop and it 
develops into shocks above z = 500 km. Figure [3] shows 
that the oscillations in the vertical velocity develop a 
clear saw-tooth shape with sudden decreases of the ve- 
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Fig. 2. — Time evolution of ^/poVl on g (left) and Y/po^trans (right) 
for the simulation with 50 s harmonic driver. Time increases from 
bottom to top, with an elapsed time between panels of a period 
(r). The image color coding is such that black represents negative 
values and white represents positive values. Horizontal lines are 
contours of constant v^/cg. 




t [min] 

Fig. 3. — Vertical velocity in the units of local sound speed at 
the axis of the sunspot at several heights in the simulations with 
50 s harmonic force. From bottom to top: z = km, z = 250 km, 
z = 500 km and z = 750 km. 

locity followed by slower increases. They present peak- 
to-peak variations of almost 8 km s _1 and their period 
is 50 s, the same period imposed by the excitation pulse. 

The second mode transformation is the acoustic fast 
mode which is transmitted as a magnetic fast mode in 
the region va > cs, where the magnetic field dominates. 




\\\m 
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RADIAL DISTANCE [Mm] 

Fig. 4. — Acoustic (top) and magnetic (bottom) flux for the simu- 
lation with a 50 s harmonic driver located at the axis of the sunspot 
averaged over the stationary stage of the simulations. Horizontal 
white line is the height where sound speed and Alfven speed are 
equal. Vertical lines are magnetic field lines. The axis are not to 
scale. 

The evolution of this magnetic mode in the first 200 s 
of the simulation is clearly seen in Fig. [2] (panels e- 
g) in the transversal velocity as the wave which moves 
away from the axis of the sunspot just above the surface 
v a ~ c 5- This mode is also visible in the transversal 
magnetic field variations (Fig. [T]i). Due to the rapid in- 
crease of the va with height the fast magnetic mode re- 
fracts and reflects back to the sub-photosphere, showing 
a b ehavior similar to the two-dim ensional case considered 
bv lKhomenko fc Colladoil (|2006h . 
When the reflected magnetic fast wave reaches again 
c| layer, it suffers two new transformations: 
a fast to slow transformation, resulting in a magnetic 
wave in the va < cs region; and a fast to fast transmis- 
sion, which produces a new acoustic wave below va = cs- 
First, we discuss the slow magnetic wave. In Fig. [T] (pan- 
els b-d) it is clearly visible in the transversal and per- 
pendicular velocities and the magnetic field variations, 
at horizontal locations inside a radius of 1.5 Mm around 
the axis and at heights between z = — 1 Mm and z = 
Mm. Observing the temporal evolution at the beginning 
of the simulation one can see that in Fig. [2^ the fast 
magnetic wave above v\ = c| has been refracted down 
and it is located between the two horizontal lines which 
indicate surfaces of constant v\/cg, but it has not arrived 
to the lower one, so the new transformation has not been 
produced. In the next time step (Fig. [2f) the wave has 
already been transformed in a slow magnetic mode in 



the v\ 



the region below v 



The wavelength of this slow 



mode decreases as the wave propagates to deeper layers 
because of the drop of the Alfven speed, which falls from 
15 km s _1 at z = Mm to 2 km s _1 at z = — 1 Mm. Due 
to the higher density at the deeper layers the amplitude 
of this wave also decreases as it propagates down. 

On the other hand, after the downward fast to fast 
transmission has occurred from the refracted fast wave, 
another fast acoustic wave appears in the region va < cs- 
It is visible in longitudinal velocity Fig. [TJi and in pres- 
sure in Fig. [T]i below the layer v\ = c|. The presence of 
this new acoustic mode can be checked comparing Fig. 
[2t and Fig. [2)3. In the latter one, there is a new wave 
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situated at both sides of the axis of the sunspot at a 
radial distance between 0.5 and 1.5 Mm which can be 
seen in longitudinal velocity. This mode appears after 
the reflected fast mode in the region va > cs reaches 
the surface v\ = c|. It propagates faster than the slow 
magnetic wave mentioned before, with a speed close to 
the sound speed, and its wavelength is larger than that 
of the slow mode. It keeps the direction of the incidence 
of the fast magnetic wave in the layer where v\ = c|, so 
it propagates down with some inclination with respect to 
the vertical, moving away from the axis. 

In order to investigate the presence of the Alfven mode 
in this simulation (either before of after the transforma- 
tions) we have plotted in Fig. the velocity component 
in the direction e per p. As expected, this velocity com- 
ponent has a node in the plane y = Mm. Thus, we 
present in Fig. a vertical cut out of this plane, at 0.4 
Mm from the center of the computational domain. In 
this projection, the Alfven wave, if present, should ap- 
pear as velocity oscillation above v\ — c| layer. The 
inspection of Fig. shows that there are no oscilla- 
tions in the magnetically dominated layers that can be 
identified an Alfven mode. Our conclusion is that when 
the wave driving occurs at the axis of the sunspot, no 
conversion to Alfven waves happens. 

3.4.3. Propagation in the upper atmosphere 

Velocity contours in Fig. [TJi show the presence of some 
waves at a radial distance between ±1.5 and ±6 Mm 
near the top of the computational domain. These waves 
are specially tricky. From the first look at the figures 
and from the time evolution of the snapshots one may 
have an impression (from the inclination of the wave 
front) that they propagate across field lines, opposite 
to the normal behavior of a slow acoustic wave. How- 
ever, we verified that their propagation speed is equal 
to cs, indicating their acoustic nature. The analysis of 
their wave numbers indicates that the direction of prop- 
agation of these waves is close to the inclination of the 
magnetic field lines. Thus, we came to the conclusion 
that these waves are slow acoustic waves created from 
the continuous transformation of the fast acoustic mode 
which moves away from the driver across the field lines 
below v\ = c| (visible in Fig. [TJd in transversal velocity 
and in Fig. [lji in pressure variations from a radial dis- 
tance of ±2 Mm to ±6 Mm) when it reaches this layer. 
As the wave front closer to the axis of the sunspot gets to 
the surface v\ = c| earlier, the slow acoustic wave that 
it produces has the wave front inclined in the direction 
to the axis. Due to this fact, it looks like its propagation 
is across the field lines. This hypothesis also explains the 
nodes located at a radial distance of ±1.5 and ±3 Mm. 

3.4.4. Acoustic and magnetic wave energy fluxes 

Figure |4] shows the acoustic and magnetic fluxes aver- 
aged over the stationary stage of the simulations. Around 
95% of the flux at the location of the driver is acoustic 
flux and it corresponds to the fast mode. At the center of 
the sunspot, where the magnetic field is almost vertical, 
the fast to slow transformation is very effective and the 
region above the layer v\ = c| at the axis is dominated 
by the slow acoustic mode. When the angle between the 
direction of propagation of the wave and the magnetic 
field is different from zero, the fast to fast transmission 
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Fig. 5. — Sound speed profile (solid line) at the axis of the 
sunspot. Diamonds represent the phase velocity of the slow acous- 
tic mode measured from the simulations with 50 s harmonic force. 

is produced and it forms the two lobes which are visible 
in the magnetic flux just above the layer v\ — c|. The 
magnetic flux reaches a height of 0.5 Mm before the fast 
waves are refracted back toward the photosphere. In the 
low-/3 region, the fast magnetic mode has an important 
contribution from z = Mm to z = 0.5 Mm for radial 
distances below 2 Mm, except at the axis of the sunspot. 
The transformation of the refracted fast magnetic wave 
when it comes back towards the photosphere generates 
acoustic as well as magnetic flux in the high-/3 region, 
corresponding to the new fast and slow modes, respec- 
tively, which propagate downwards. As expected, we do 
not find any propagating magnetic flux along the field 
lines, that may be associated with an Alfven mode. 

3.4.5. Slow acoustic mode phase velocities 

The phase velocity of a linear high-frequency slow 
mode wave in a magnetically dominated region is equal to 
cs- In our simulation, a slow acoustic wave appears above 
the v\ = c| layer produced after the mode transforma- 
tion. This wave allows us to check whether or not the 
velocity of waves involved in the simulations in this com- 
plex sunspot model corresponds to that expected from 
theoretical considerations. 

Figure [5] presents the results of this test. The solid 
line in Figure [5] shows the stratification of the sound 
speed with height at the sunspot axis and diamonds in- 
dicate the phase velocity of the slow mode wave mea- 
sured at each grid point from the simulations. Note that 
diamonds are only plotted at heights above z = —200 
km, after the mode transformation has been completed. 
The velocity of the slow wave matches well the local 
sound speed at heights from z = —200 km to z = 700 
km. Higher than z = 700 km the wave starts to prop- 
agate faster than cs, since the velocity amplitude of 
the wave approaches the sound speed and non-linearities 
start playing an important role. 

3.4.6. Slow acoustic mode amplitudes 
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Fig. 6. — Amplitude of the slow acoustic mode vertical velocity 
(solid line) in the simulations with 50 s harmonic force. Dashed 
lines gives t he analytical curve for an acoustic-gravity wave with 
50 s period (Mihalas &; Mihalas 1983). 

As the slow acoustic wave propagates up in the atmo- 
sphere of the sunspot the amplitude of the velocity oscil- 
lations increases due to the density drop. The kinetic en- 
ergy of this wave is proportional to pv 2 and must be con- 
served. Thus, a decrease of the density must be accom- 
panied with an increase of the velocity. The atmosphere 
of the sunspot is very complex and it includes vertical 
and horizontal gradients in all the magnitudes. Because 
of that there is no analytical expression for the variation 
of the amplitudes of the waves with height. However, 
we can compare the particular case of the wavefront of 
the slow acoustic mode wave which propagates along the 
axis of the sunspot with a case of a linear acoustic wave 
which propagates upwards in a gravitationally stratified 
atmosphere permeated by a magnetic field parallel to the 
direction of gravity. In this theoretical case, the ampli- 
tude of the wave is given by 



A(z) = A exp^j —J, (37) 

In Fig. [6]we show the amplitude of the vertical velocity 
at each grid point at the axis of the sunspot (solid line) 
for the stationary stage of the simulations. The force is 
located at z = —500 km and the amplitude of the wave 
decreases until it reaches z = —200 km. This initial de- 
crease is due to the part of the energy of the source that 
goes into other wave types because of the mode trans- 
formations. More or less at this layer the wave is trans- 
formed into a slow acoustic mode as it propagates up 
while its amplitude increases. When the wave reaches 
the height z = 850 km its amplitude drops very fast as 
a consequence of the large diffusivity that was imposed 
at high layers in order to stabilize the numerical sim- 
ulation. In this figure is also overplotted the expected 
amplitude according to Eq. (|37|) (dashed line), starting 
from the height where the wave has already been trans- 
formed into a slow acoustic mode in the region va > cs- 



From z = —200 km to z = 700 km the numerical ampli- 
tude agrees with the analytical one, while from z = 700 
km to z = 850 km the numerical amplitude is lower than 
the analytical one. This happens because the wave de- 
velops into weak shocks and the linear approximation 
for the amplitude increase is no longer valid. Note that 
the amplitudes of the oscillations (Fig. [6]) as well as the 
phase velocity (Fig. [5j) show discrepancies with the linear 
theory at the same heights. 

3.5. Case of 50 s harmonic force located off the axis 

Fig. [71 gives a three-dimensional view of the vertical 
velocity in the simulation run with 50 s harmonic force 
located at x = — 3 Mm off the sunspot axis. This figure 
clearly shows the asymmetry of the wave front with re- 
spect to the axis. In the lower part of the domain, the 
fast (acoustic) waves can be appreciated propagating in 
circles away from the source with a visibly lower ampli- 
tude toward the axis. In the upper part of the domain, 
slow (acoustic) waves are the dominating ones, propagat- 
ing along the inclined magnetic field lines. 

Fig. [5] presents the projected velocities in the three 
characteristic directions at the horizontal cut of the 
simulation domain taken at the middle photosphere at 
z = 300 km. Fig. [9]shows the snapshots of some variables 
in the vertical cut through the domain. Both figures cor- 
respond to the same time moment of the simulations at 
t = 820 s. These simulations have many features in com- 
mon with the previously considered case of the driving 
force located at the sunspot axis. A set of fast (acous- 
tic) and slow (magnetic) modes is generated below the 
layer v\ = c|, propagating upwards and suffering sev- 
eral transformation after reaching this height. Similar 
to the previous case, slow (acoustic) and fast (magnetic) 
modes are produced after the mode transformation in the 
magnetically dominated upper atmosphere. The conver- 
sion to slow and fast modes in the low-/? region only 
presents slight changes in comparison with the simula- 
tion with the driver placed at the axis. One of these 
changes is the presence of an asymmetry with respect to 
the axis. For example, the transformation of the down- 
ward propagating refracted fast (magnetic) mode into 
the slow (magnetic) mode below the surface (Fig. [9j pan- 
els b-d) presents such asymmetry. Due to the particular 
combination of the field inclination and the direction of 
propagation of the refracted fast mode, this transforma- 
tion is much more efficient on the right from the source in 
the direction toward the axis (locations between x = — 3 
and -1 Mm at z « -0.5 Mm). 

The most important changes are present in the velocity 
component in the direction e perp (Fig. |5t). At height 
of z = 300 km, this component shows variations that 
do not correspond to the fast (magnetic) mode. The 
latter has been already reflected down. As demonstrated 
by the analysis of the magnetic energy flux below, the 
variations observed in the snapshot of e pe r P at heights 
above z ~ 300 km at locations between x = — 4 -i — 2 
Mm and y = —1 + 1 Mm correspond to the Alfven wave. 

The acoustic and magnetic fluxes are shown in Fig. 
[TOl These fluxes, in general, show a pattern similar to 
the case of driving at the axis, but some important asym- 
metry is also present. The acoustic flux of the slow mode 
in the low-/? region is oriented in the direction of the field 
lines and it is slightly lower than in the previous simu- 
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Fig. 7. — Variations of yfpv z at an elapsed time t=820 s after the beginning of the simulations for the 50 s harmonic force at 3 Mm from 
the axis of the sunspot. Grey inclined lines are magnetic field lines. Blue colors represent upward movement while orange/red colors are 
downward movement. 
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Fig. 8. — Variations of the velocity in the direction ei ong (a), etrans (b), and e per p (c), all of them scaled with factor ^/po of the unperturbed 
density, at z = 300 km and at an elapsed time t = 820 s after the beginning of the simulations for the 50 s harmonic force located at x = — 3 
Mm , y = Mm and z = —0.5 Mm. The image color coding is such that blue colors represent lower values and red colors are higher values 
with respect to the mean. All images have the same scale. Concentric dashed lines are contours of equal magnetic field. 



lation due to the larger angle between the direction of 
propagation of the fast mode before reaching the layer 
v a = c l an d the magnetic field ((f). An important frac- 
tion of the magnetic flux appears in the lobe located 
at v\ = c| line in the direction toward the axis of the 
sunspot present due to more effic ient fast to fast mode 
transmission with increasing cp ([Callyl [2QQ5h . Compared 
to Fig. H] there is also more magnetic flux present in the 
upper part of the atmosphere above v\ = c|, apparently 
directed along the magnetic field lines. It is not clear, 
however, whether this flux corresponds to fast, not yet 
completely reflected wave, or to the Alfven wave. To 
clarify this issue we have calculated the magnetic flux 
following Eq. ([36]) . but using only velocity and magnetic 
field projections in the direction e per p. We expect that 
the longitudinal projection of this quantity gives us in- 
dications about the presence of the propagating Alfven 
waves. 

Figure fTTI illustrates the result. Top panel corresponds 
to a vertical cut in the plane y = —0.4 Mm, normalized 
at every height to its maximum value at this height. The 
bottom panel is a horizontal cut in the plane z = 0.9 Mm. 
The white colors (positive flux) mean upward energy 



propagation, while the black colors (negative flux) mean 
downward energy propagation. Indeed these plots reveal 
the energy flux associated to the Alfven wave, which 
clearly propagates upwards along the field lines. The 
Alfven mode has a node at the plane y = Mm, where 
the driver was located, so conversion to this mode is only 
produced when the wave vector forms a certain angle 
with the magnetic field (different from zero). This result 
is in qualitative agreement with the r ecent investigation 
of th e conversion to Alfven waves by (jCallv fc Goossensl 
2008). However, even at the location where the contri- 
bution of the Alfven wave energy flux to the total energy 
flux is maximum, its flux is still around 20 times lower 
than the acoustic flux at this location. So it means that 
in the sunspot magnetic field configuration and for the 
driver location considered here, the transformation from 
fast (acoustic) to Alfven wave is much less effective than 
the transformation to the slow (acoustic) wave in the 
magnetically dominated upper atmosphere. 

3.6. Case of 180 s harmonic force located at the axis 

This simulation shows many features which also ap- 
peared in the 50 s harmonic case, since the response of 
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Fig. 9. — Variations of the velocity in the direction ei ong (a), etrans (b), and e per p (c), all of them scaled with factor ^/po of the unperturbed 
density; magnetic field in the direction etrans, (d) and pressure (e) at an elapsed time t=820 s after the beginning of the simulations for 
the 50 s harmonic force located at x = —3 Mm off the sunspot axis. The format of the figure is the same as Fig. ^ 

the magnetic atmosphere is quite similar at both frequen- 
cies because they are above the cut-off frequency. So, in 
this case we omit the discussion of the mode transforma- 
tion and only show results of the acoustic and magnetic 
flux calculations. These fluxes are presented in Fig. [121 
This figure is very similar to Fig. 2J except for much 
larger wavelength of the fluctuations, in agreement with 
the larger temporal period of waves. The average acous- 
tic flux shows that at the axis of the sunspot the fast to 
slow mode transformation is also effective for waves with 
180 s period. After the transformation, slow (acoustic) 
waves propagate acoustic energy upwards. The fast mode 
visible in the magnetic flux above the v\ = c| also has 
longer wavelenghts than in the 50 s harmonic simulation. 
The magnetic flux vanishes at the axis of the sunspot at 
high layers. Due to the reflection of the fast (magnetic) 
wave, there is no magnetic flux above the certain height, 
also away from the axis. In this simulation we find no 
traces of the Alfven mode. 





RADIAL DISTANCE [Mm] 

Fig. 10. — Average acoustic (top) and magnetic (bottom) flux for 
the 3D simulation with a 50 s harmonic driver located at 3 Mm 
from the axis of the sunspot averaged over the stationary stage of 
the simulations. The format is the same as Fig. 0] 



3.7. Case of 300 s wavelet force located at the axis 
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Fig. 11. — Magnetic flux of the Alfven mode. Top: Vertical cut 
in the plane y = —0.4 Mm, normalized at every height. Vertical 
white lines are magnetic field lines and horizontal white line is the 
layer where Cg = v\- Bottom: Horizontal cut in the plane z = 0.9 
Mm. Thin dashed lines are contours of equal magnetic field. In 
both panels, thick dashed lines mark the location of the other plot. 



Figure [13] shows vertical snapshots of several variables 
after 820 s of the simulation with the driver emitting a 
spectrum of waves with a central frequency at 3.33 mHz 
given by the Eq. (|29|) . 

According to the stratification of the atmosphere, at 
the axis of the sunspot at z = —700 km the cut-off fre- 
quency is v c = 3.3 mHz. It increases with height reaching 
a maximum at the height of z = km, where its value is 
v c = 5.8 mHz. It means that 5 minute acoustic waves can 
not propagate upwards above z = Mm, since they are 
evanescent in the vertical direction, and can only prop- 
agate horizontally. Therefore, the cut-off frequency is a 
critical value for the wave propagation in this case, and 
the behavior of waves below and above the cut-off layer 
is expected to be completely different. 

The source located at z = —500 km drives waves with 
frequencies below as well as above z/ c , but most of its en- 
ergy goes to the band around 3.33 mHz. This generates 
a fast acoustic wave with an amplitude below 100 m s _1 
and 5 minute period, which can propagate upwards only 
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Fig. 12. — Acoustic (top) and magnetic (bottom) flux for the 
simulation with a 3 min harmonic driver located at the axis of the 
sunspot averaged over the stationary stage of the simulations. The 
format is the same as Fig. [4] 

until the height z = Mm. Waves with this frequency 
are evanescent at higher layers and their vertical wave- 
length occupies the whole upper part of the simulated 
atmosphere. The amplitude of their longitudinal veloc- 
ity slightly increases with height. 

3.7.1. Three-dimensional mode transformation 

Part of the energy of the driver which reaches the sur- 
face v\ = c| is transformed into a fast wave above this 
height. Due to its magnetic nature it becomes unaffected 
by the cut-off frequency. As in the previous simulations, 
the transversal velocity in Fig. [T3b and the transversal 
magnetic field in Fig. [Kill show that the fast magnetic 
wave in the region va > cs is reflected because of the gra- 
dients in the Alfven speed. Once the wave comes back to 
the sub-photospheric layers below v\ = c|, it keeps its 
magnetic nature and propagates downwards in a form of 
a slow wave with decreasing wavelength due to the drop 
of the Alfven speed, but keeping its 5 minute period. 

The maximum cut-off frequency at the axis in this 
sunspot model is v c = 5.8 mHz, so waves with higher 
frequencies can still propagate upwards through the at- 
mosphere. The fast acoustic modes generated by the 
driver with frequencies higher than v c are transformed 
into propagating slow acoustic modes in the region above 
v a ~ c 5- The contours in Fig. [T3k represent constant 
longitudinal velocity. At a height around z = 900 km, the 
longitudinal velocity has maximum power at 3 minute 
period, which corresponds to the frequency above v c re- 
ceiving more energy from the driver. 

3.7.2. Frequency change with height 

The increase of the amplitude of 3-minute waves ac- 
cording to the drop of the density (compared to much 
weaker increase of the amplitude of the evanescent 5- 
minute waves) leads to the power spectrum at chromo- 
spheric heights dominated by 3-minute waves. There, 
their amplitudes reach almost 400 ms _1 . They do not 
develop into a saw-tooth waves because at the photo- 
sphere the driver generates low power at this frequency 
band and their amplitude increase is not enough to pro- 
duce significant non-linearities. 
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Fig. 13. — Variations of the velocity in the direction ei ong (a), etrans (b), and e per p (c), all of them scaled with factor yfpo of the 
unperturbed density; magnetic field in the direction etrans (d), and pressure (e) at an elapsed time t=820 s after the beginning of the 
simulations for the 300 s wavelet force located at x = km, y = km and z = —500 km. The format of the figure is the same as Fig. ^ 
Violet lines represent contours of equal cut-off frequency. The inner one is v c = 5.6 mHz and the outer one is v c = 4 mHz. 



Figure [TH shows the power spectra at different heights, 
from z = —250 to z = 750 km. At z = 750 km it shows 
a clear peak at v = 5.9 mHz, corresponding to a period 
of 170 s. At the lower layers, a peak below 4 mHz dom- 
inates. One can see in this figure how the amplitude of 
the low- frequency peak increases with height. However, 
this increase is weaker than the one of the peak at high- 
frequency around v « 6 mHz. Due to this behaviour, 
the high layers are dominated by oscillations with period 
around 3 minutes. 

A simulation with a 300 s harmonic driver was also 
performed. In this case high layers develop evanescent 
5 minutes period waves, but there is no trace of 3 min- 
utes oscillations at the chromosphere, so power in this 
frequency band can not be produced if these frequencies 
are not excited by the driver. Thus, we can conclude 
that the mechanism that produces the change in fre- 
quency of oscillations in the umbra from the photosphere 
to the chromosphere is the linear propagation of waves 
with 3 minute power which come directly from the pho- 
tosphere and dominate over the evanescent long period 



waves. This conclusion goes in line with the results of the 
observational study of sunspot waves simu ltaneously at 
the ph otosphere and the chromosphere bv lCenteno et al.l 
([2006h . 

3.7.3. Acoustic and magnetic wave energy fluxes 

The magnitudes ^/pov\ ong , y/pov tran and y/pov perp plot- 
ted in Fig. [T3l (panels a-c), respectively, show that most 
of the kinetic energy remains in the photosphere and be- 
low. Most of the energy introduced by the driver goes 
into the waves in the 5 minute band. Propagating to 
higher layers, they form evanescent waves or are trans- 
formed into fast magnetic mode waves. The first ones 
do not carry energy, while the second ones are reflected 
back to the photosphere. Thus, waves in the 5 minute 
band can not supply energy to the chromosphere, if the 
driving force is located at the sunspot center. 

This case differs from the simulations with shorter pe- 
riods described in Sects. 13. 4| I3.5l and l3.6| where the slow 
acoustic wave transports to the high layers part of the en- 
ergy injected by the driver (or, in the case of 50 s off-axis 
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Fig. 14. — Power spectra at different heights at the axis of the 
sunspot for the simulations with 300 s wavelet force, normalized to 
the maximum power at the highest height. From bottom to top: 
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Fig. 15. — Acoustic (top) and magnetic (bottom) flux for the 
simulation with a 5 min wavelet driver located at the axis of the 
sunspot averaged over the stationary stage of the simulations. The 
format is the same as Fig. 3] Violet lines represent contours of 
equal cut-off frequency. The inner one is v c = 5.6 mHz and the 
outer one is v c = 4 mHz. 

driver, a smaller part of the energy is also transported 
upwards in the form of the Alfven wave). In the sim- 
ulation with the 5 minute wavelet force located at the 
axis, only the waves with frequencies higher than v c can 
provide energy to the chromosphere, and they represent 
a small fraction of the energy introduced by the driver. 

The acoustic flux in Fig. [151 shows that most of the en- 
ergy keeps below the layer c| = v\. The wavelet mainly 
drives 5 min power in a fast acoustic mode, and since 
it is evanescent, it does not propagate energy upwards 
and this 5 min power is distributed horizontally. Only 
waves with frequency higher than the cut-off frequency 
are transformed into slow acoustic modes in the low-/3 
region and carry energy to the chromosphere. They cor- 
respond to the low acoustic flux which appears at the 



center of the sunspot between z = Mm and z = 1 
Mm. The high-/3 region contains the energy of the slow 
magnetic modes generated from the secondary transfor- 
mation of the reflected fast magnetic modes. In the mag- 
netic flux, the energy of these slow modes appears in red 
for radial distances below 6 Mm. Note that, interest- 
ingly, the maximum of this energy is not located just at 
the c| = v\ line it was in the case of 50 s simulations, 
but it located below this line. This new location follows 
the line of constant v c = 4 mHz meaning that the cut-off 
effects influence also the penetration of the magnetic en- 
ergy into the higher layers. The fast mode waves in the 
high-/3 region contribute to the high acoustic flux there. 
Like in all the simulations with the driver located at the 
center of the sunspot, magnetic flux along field lines is 
negligible, and there is no conversion to Alfven waves. 



3.7.4. Propagation in the upper atmosphere 

Contours in Fig. fT3li in the region above v\ = c| (up- 
per part of the domain) far from the axis of the sunspot 
show longitudinal waves with 5 minute period which ap- 
parently mo ve ac ross field lines at the sound speed. Sim- 
ilar to Sect. 13. 4| we tried to analyze their wave number 
and phase velocity behavior with height. However, the 
wavelength of these waves in the upper layers is compa- 
rable to the size of our simulation domain in the vertical 



direction (r 



1 Mm above v\ 



c 2 s ). Because of that, 



we can not be completely sure if these waves propagate 
along or across the field lines. 

Another difficulty to understand the behavior of these 
waves lies in the fact that their frequency is below the 
cut-off frequency v c of the atmosphere. In principle, the 
cut-off frequency in the magnetically dominated atmo- 
sphere is lowered for the acoustic wave propagating along 
inclined field lines. We calculated the cut-off frequency 
taking into account the effect of the field inclination. We 
found this inclination to be insufficient to reduce the cut- 
off frequency enough to allow for the propagation of the 
5-minute waves. Simulations in a larger spatial domain 
(both in the horizontal and vertical direction) will be 
needed in the future to clarify the nature of these waves. 

4. CONCLUSIONS 

In this paper, we have addressed the problem of the 
three-dimensional wave propagation and mode transfor- 
mation of the MHD waves in the upper atmosphere 
(photosphere and chromosphere) of a sunspot model, by 
means of numerical simulations. We have presented our 
code for the calculation of the response of a magneto- 
static structure in equilibrium to an arbitrary perturba- 
tion. This code is specially designed to model the wave 
propagation in solar magnetic structures with non-trivial 
magnetic field configurations. We performed several tests 
proving the correct code performance. The comparison 
of the numerical tests with the analytical solutions (when 
they are available) or with the results from other codes 
(in other cases) shows a precise agreement and demon- 
strates the robustness of the numerical method applied 
in this code according to several aspects: the velocities of 
wave propagation, the resolution of shocks and the con- 
servation of energy. The performance of the artificial dif- 
fusivity makes it possible to resolve strong discontinuities 
with a few grid points and without producing damping of 
the solution, and the energy is well conserved. The PML 
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boundary condition allows us to calculate long tempo- 
ral series of simulations without reflections of the waves 
reaching the top boundary. 

We have presented the analysis of several simulations 
where the sunspot atmosphere was perturbed by dif- 
ferent pulses, varying their location and temporal be- 
haviour. The simulations of short period waves in the 
three-dimensional sunspot model clearly show several 
phenomena that are predicted by wave theory. We con- 
firmed that our code correctly describes the propagation 
of slow and fast modes, both in regions dominated by 
the magnetic field or the gas pressure. Our main find- 
ings, learned from the simulations, can be summarized 
as follows: 

• The conversion between fast and slow magneto- 
acoustic waves happens in three dimensions in a 
qualitatively similar way as in two dimensions. 
Waves with frequencies down to the cut-off fre- 
quency behave in the same way. The driver lo- 
cated in the gas pressure dominated region gener- 
ates mostly the fast (acoustic) mode. This mode, 
propagating to the upper layers, is transformed 
at the height where cs = va- After the trans- 
formation, a slow acoustic mode propagates up- 
wards along the field lines in the magnetically dom- 
inated atmosphere. The fast magnetic mode under- 
goes refraction and it is reflected back to the sub- 
photosphere. When it reaches again the surface 
cs — va, new transformations take place produc- 
ing another fast acoustic and slow magnetic modes 
in the region va < cs- 

• High-frequency field-aligned propagating acoustic 
waves are constantly produced in the upper mag- 
netically dominated atmosphere at locations away 
from the source. These waves appear due to the 
continuous transformation from the fast (acoustic) 
waves moving horizontally, across the field lines, 
away from the source in the gas pressure dominated 
region. On their way, the fast waves constantly 
touch the cs = va layer producing slow (acous- 
tic) waves in the upper atmosphere in the horizon- 
tal locations far from the source. We observe this 
behaviour in all simulations with different driving 
frequencies and source position. 

• The 3D simulations allow us to identify an Alfven 
mode. This mode appears only in the simulation 
with the source located away from the sunspot axis. 
It is produced after the transformation from the 
fast (acoustic) mode. We find that the transforma- 
tion efficiency from the fast to the Alfven mode is 
much lower than that from the fast to the slow 
mode. The eventual energy of the Alfven wave 
in the magnetically dominated region is 20 times 
lower than that of the slow mode. In the simu- 
lations with the driver located at the axis of the 
sunspot, where the angle between the direction of 
the upwards propagating wave and the magnetic 
field is zero, we find no indications of the transfor- 
mation to the Alfven mode. 

• The analysis of the wave energy fluxes suggests 
that in the high-frequency cases (above the cut-off) 



the wave energy can reach the upper atmosphere 
most efficiently in the form of slow (acoustic) field 
aligned propagating waves. After some height in 
the middle photosphere there is no magnetic flux 
corresponding to the fast (magnetic) waves as their 
energy is reflected. If the driver is located away 
from the axis, some small part of the energy also 
can propagate upwards in the form of an Alfven 
wave. 

• Both magnetic and acoustic energy of the low- 
frequency waves (smaller than the cut-off) remain 
almost completely below the level cs = va and 
do not reach upper layers. This happens because 
the energy of the fast magnetic modes in the up- 
per layers is reflected, and the evanescent acoustic 
slow modes do not propagate any energy at these 
frequencies. The comparison between the magnetic 
fluxes at high and low frequencies shows that the 
magnetic flux reaches smaller heights for the low- 
frequency waves. 

• When the driver excites a spectrum of waves, we 
observe a change of the dominant frequency of os- 
cillations with height from the photosphere to the 
chromosphere. The driver excites a spectrum with 
a maximum power at 3.33 mHz frequency. Waves 
at this frequency are evanescent in the atmosphere. 
At higher frequencies (above v c — 5.8 mHz in our 
sunspot model) the waves can propagate upwards 
along the field lines. Due to the larger amplitude 
increase with height of the propagating waves, com- 
pared to the evanescent waves, the 3-minute waves 
(y « 5.9 mHz) dominate the power spectrum in 
the chromosphere. This behavior, obtained in the 
simulations, is similar to the observed one. 

One of the questions that arises from these results 
is the evaluation of the validity of the decomposition 
performed to separate the fast mode from the Alfven 
mode in the low-b e ta pla sma. The decoupling following 
iCallv fc Goossensl (|2QQ8[ ) is valid for an idealized case 
when considering a uniform magnetic field, for a plane 
wave with a constant wavenumber perpendicular to grav- 
ity, and it is obtained asymptotically in the limit of in- 
finite Alfven speed. Although the realistic atmosphere 
used in our calculations does not fulfill these restrictions, 
a coherent picture is retrieved from the magnetic flux of 
the velocity and magnetic field components in the direc- 
tion given by Eq. [3H showing upwards propagation along 
magnetic field lines. The result that naturally emerges 
from this decomposition justifies the application of the 
method. 

Another issue that complicates the decoupling of both 
modes is the excitation of small horizontal wavenumbers 
due to the limited horizontal extent of the driver. The 
fast mode is refracted due to the rapid increase of the 
Alfven speed with height, but the altitude at which it 
happens depends on the wavenumber, being higher for 
lower wavenumbers. Since most of the power excited 
by our driver lies in the range of wavenumbers below 
~ 1/R scr: some of the fast mode waves might still be 
partially refracted in the limited 1 Mm atmosphere above 
the cs = va surface. The magnetic flux of these waves 
may reach high layers, complicating the separation of 
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the fast and Alfven modes. In this scenario, we would 
expect a continuous transition between the upward lon- 
gitudinal magnetic flux and the refracted one at different 
heights where all these waves are being refracted. How- 
ever, the magnetic flux of the Alfven wave that we re- 
trieve is clearly delimited along magnetic field lines, con- 
firming that this flux mostly corresponds to the Alfven 
mode. 

The most important achievement reached by the 
development of our numerical code is the possibility to 
investigate the three-dimensional mode transformation 
in realistic conditions imitating a sunspot atmosphere. 
This, together with the possibility to study large-period 
waves in the layers where they are observed, gives an 
opportunity for the direct comparison between our 
numerical simulations and solar spectropolarimetric 
observations. Simulations of sunspot high layers repre- 
sent a hard challenge due to the exponential increase 
of the Alfven speed with height. Our sunspot model 
presents an Alfven speed of almost 1000 km s _1 close 
to the upper boundary of the domain, limiting the 
time step and making the calculations very expen- 
sive. Despite this, our code manages to describe the 
waves well, including the correct performance of the 
boundary PML layer. Note that other works on waves 
in non-trivial magnetic configurations have been re- 
stricted either to two-dimensional h i gh-frequency cases 
(iCargill et~aD 119971: [Rosentha l et al.l l2QQ2t iHasan etaLl 
12003 iBogdan et al.ll200l iKhomenko fc Colladosl I2QQ6D 
or to the study of helioseismic waves, w here the problem 

of high Alfven speed is avoid ed (Cally & B ogdanl 

19971: [Parchevs ky fe Kosovicheyl 120091; iHanasogel 
2008; Cam eron et al.l 120081: iMoradi et al.l 120091* 



Khomenko et aLl l2QQ9). 



The strategy applied in our code allows the direct com- 
parison with observations by means of spectral synthe- 
sis. The simulations presented in the paper reproduce 
the frequency change with height, indeed o bserved in the 
sunspot atmospheres (Ce nteno et a l. 2006). In the future 
we plan to perform a more detailed comparison with solar 
data. This will be done by exciting the sunspot magneto- 
static model in equilibrium with velocities obtained from 



spectropolarimetric observations, and by comparison of 
the simulated wave parameters in the photosphere and 
chromosphere with those obtained from simultaneous ob- 
servations in different spectral lines. 

Perhaps the most interesting simulation considered in 
our paper is the one with the source exciting a spec- 
trum of waves close to the solar one. This case has 
a special relevance because theoretical models of wave 
transformation, existing as of today, are best valid in 
the high-frequency limit and do not address the behav- 
ior of waves at frequenc ies below the cut-off frequency 
(jSchunker fc Callvll2006f ). Our initial results show that 
almost no energy of the 5-minute waves propagates into 
the higher layers, at least in the situation of the sunspot 
model and source located at the axis considered in the pa- 
per. These results need further exploration. In our future 
work we will perform simulations of 5-minute waves in 
larger computational boxes, varying the excitation mech- 
anisms (location and properties of the source), as well as 
the sunspot model. In particular, an interesting question 
is under which conditions 5-minute Alfven waves can be 
excited by the mode transformation. These waves can 
still propagate some energy into the upper atmosphere 
of solar active regions, and thus understanding the con- 
ditions of transformation to these waves and their ener- 
getics in sunspots is important. The number of works in 
the literature with numerical calculations including the 
transformation to Alfven waves is scarce. The mos t rel- 
evant study is the one by Cal ly fc Goossensl ((2008) who 
find that the transformation to an Alfven mode is effec- 
tive at certain angles of inclination and azimuth of the 
magnetic field. 
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through projects AYA2007-63881 and AYA2007-66502. 
The simulations have been done on LaPalma supercom- 
puter at Centro de Astrofisica de La Palma and on 
MareNostrum supercomputer at Barcelona Supercom- 
puting Center (the nodes of Spanish National Supercom- 
puting Center). To represent the data we have used VA- 
POR software ( |http://www.vapor. ucar.edu] ). 
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Fig. 16. — Standard Riemann shock tube problem at t = 0.2. Crosses represent numerical solution, while lines represent the analytic 
solution with an exact nonlinear Riemann solver. Inner plot on left bottom panel helps to visualize better the discontinuity at x = 0.68. 
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APPENDIX 
TESTS ON NUMERICAL PERFORMANCE 
In this appendix we describe the results of standard numerical tests to verify the code performance. 

ID Riemann shock tube test 

The Riemann shock tube test (Sod 1978; Caunt fc Kor pi 2001) has been simulated in order to test the behavior of 
the hydrodynamical part of the code, including discontinuities in the properties of the fluid. This allow us to analyze 
how the artificial viscosity copes with shock capturing. The physical domain size is unity and the initial conditions 
include a discontinuity at x = 0.5. On the left we have density p\ = 1 and pressure pi = 1, while on the right 
pi = 0.125 and p\ = 0.1. The ratio of specific heats is 7 = 1.4 and the initial velocity as well as the magnetic field 
are set to zero. The problem has been simulated in one dimension with a resolution of 256 grid points and closed 
boundaries. 

Figure [16] shows the density, velocity, internal energy density per unit mass and pressure at time t=0.2 for the 
simulation compared to the analytical solution. From left to right the plot shows a rarefaction wave (from x = 0.25 
to x = 0.5), a contact discontinuity (at x = 0.68) and a shock front (x = 0.85). The position of all of them matches 
precisely with the analytical solution and the magnitudes of the fluid properties are correct. The contact discontinuity 
for the energy was inevitably smoothed, as shown in the plot at the bottom left where a region around it is amplified, 
but still a 93% of the amplitude of the discontinuity is covered with 10 grid points. Moreover, the shock front is 
resolved with 3 grid points, which proves the good performance of the code in shock capturing. 




1.5D Brio & Wu shock tube 

To test th e formation of mag net ohydro dynamic shock waves we use the MH D analog of the Sod shock tube problem 
descr i bed bvlBrio fc Wul (1 19881) . which has been widely used in previous works (jStone fc Norm an 1992; C aunt fc Korpil 
120011 : IShelvag et al.ll2008h . We can compare our results with those given in the literature, since no known analytical 
solution exists for the evolution of this problem. In this 2D test the fluid is initialized in a physical domain from x = 
to x = 1 and with a discontinuity in density, pressure and magnetic field normal to the direction of motion located at 
x = 0.5. Parameters at the left hand side from the discontinuity are p\ = 1, p\ = 1 and B y \ = yT^o? an d at the right 
hand side are p 2 = 0.1, p 2 = 0.1 and B y2 = — ^/P0' All the domain is permeated with a constant magnetic field along 
the direction of motion B x = 0.75y7^o and the adiabatic index 7 is set to 2. In this case the resolution is 800 grid 
points in the direction of the shock wave propagation, similar to the other published works. All the boundaries are 
closed. 

The density, velocity in x direction, pressure and magnetic field in y direction are shown in Figure [TTJ This MHD 
Riemann problem produces a complex solution with several components: the waves moving to the left are a fast 
rarefaction wave and a slow compound wave (consisting of a slow rarefaction attached to a slow shock), and the waves 
moving to the right include a contact discontinuity, a slow shock and a fast rarefaction wave. Comparison of our 
results with other works show that the waves have propagated with the correct velocities and have similar magnitudes, 
indicating that they are in good agreement with other solutions. The slow shock is resolved again with only 3 grid 
points. 

2D Orszdg-Tang vortex 

The next test is the O rszag- Tang vortex, which was originally s t udied by Orszag & Tang (1979) and has been use d 
to probe several codes (jRvu et al.l 119951 : iDai fc Woodwardl 119981 : lLondrillo fc Del Zannal l2000: She lvag et aI1l2008[ ). 
This problem allows us to demonstrate the robustness of the numerical scheme used in our code solving the two- 
dimensional interaction of non-linear shock- waves and also to compare qualitatively the code with other codes. The 
initial conditions for density and gas pressure are constant, with p = 25/(367r) and p = 5/(127r), the magnetic field 
B x = — sin(27n/) and B y = sin(47nr) and the initial velocity v x = — sin(27n/) and v y = sin(27nr). Therefore, the initial 
flow is a velocity vortex superimposed to a magnetic vortex, with a common X-point, but with different structure. The 
initial Mach number is Mo = 1 and the adiabatic index is set to 7 = 5/3. In our simulation for this problem we have 
chosen a unit size in horizontal and vertical dimensions and the resolution of the simulation box is set to 512 x 512 
grid points. Figure [18] presents the density at time t = 0.5, showing precise agreement with the other published works. 



3D Acoustic wave 

Since our code is oriented for the simulations of waves, it is necessary to test how well it can approximate the known 
analytical solutions for different types of waves in a stratified atmosphere, their propagation speeds, amplitudes and 
shock development. The analytical solution of an acoustic wave propagating in an isothermal at mosphere with vertical 
strati fication due to gravity and permeated with a constant vertical magnetic field is known from iFerraro fc Plumptonl 
(1958). The vertical velocity for wave with frequencies uj above the cut-off frequency uj c = / yg/(2cs) follows as: 



v z (z, t) = De 



z/2H , 



cs 



-z + cut 



(Al) 
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Fig. 18. — Mass density at t=0.5 of the Orszag-Tang vortex simulation. 

where Ho = c|/ (7*7) is the pressure scale height and cs is the sound speed. This analytical solution with a period of 
15 s and starting amplitude of 10 m s _1 has been introduced as a bottom boundary condition without variations in 
the horizontal directions and its evolution in time has been calculated with the numerical code. The physical domain 
is set to 1000 km in the vertical direction with a resolution of 200 grid points, using 20 of them as PML domain at the 
top. In both horizontal dimensions the computational domain consists of 7 grid points, covering 35 km, and periodic 
boundary conditions were used in these directions. In Fig. [19] we show the comparison between the numerical and the 
analytical linear solution for the vertical velocity after 197 s of simulations. The numerical solution matches the exact 
solution in the computational domain, while it is damped effectively by the PML layer. It is important to note that 
the amplitude and propagation speeds are both described correctly by the numerical solution and that the numerical 
diffusion does not affect the amplitude increase with height. No spurious reflections are present showing the correct 
action of the PML boundary in this example. 

The dashed line represents the difference between both solutions, which also increases with height. The main 
contribution to this difference are the nonlinear terms that are taken into account in the numerical calculation, but not 
in the analytical solution. We have checked that reducing the driver amplitude by a factor of two reduces the difference 
by exactly a factor of four. The wavelength of the difference is twice shorter compared to that of the wave. These two 
arguments prove that the difference between the analytical and numerical solutions is mostly due to non-linear terms. 

3D Alfven wave 

As a next step, the response of the numerical scheme to the propagation of an Alfven wave in an isothermal, stratified 
atmosphere with vertical magnetic field is analyz ed. The analytical solution was developed by iFerraro fc Plumptonl 
(|1958l ). and, according to iKhomenko et aTl (j2QQ3h . the solution for the horizontal velocity can be written as 



where Jo and Yq are Bessel functions and va is the Alfven speed. The same atmosphere as in the previous section was 
used in this test, but now the bottom layers were excited with the solution of an Alfven wave of period 10 s and an 
amplitude of 10 ms" 1 as a boundary condition. This driver generates the propagation of Alfven waves toward higher 
layers of the atmosphere. The horizontal velocity of both the numerical and exact solutions for time t = 115 s is shown 
in Fig. [20j demonstrating a very good match. 

The dashed line shows the difference between the numerical and the analytical solution. In this case, the nonlinearities 
are not so important since the amplitude is lower than the one of the acoustic wave, and the discrepancy between 
both solutions is due to the reflection produced at the top boundary. The PML layer results more problematic for 
transversal waves which oscillate parallel to the interface between the PML media and the physical domain, and they 
give rise to reflections of 5-6% of the velocity. 
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Fig. 19. — Vertical velocity profile of vertically propagating acoustic waves in an isothermal, stratified atmosphere with vertical constant 
magnetic field at t=197 s. Solid line: exact solution; diamonds: numerical solution. Dashed line is the difference between both solutions. 
The vertical dashed line indicates the location of the PML interface. 
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Fig. 20. — Horizontal velocity profile of Alfven waves in an isothermal, stratified atmosphere with vertical constant magnetic field at 
t=115 s. Solid line: exact solution; diamonds: numerical solution. Dashed line is the difference between both solutions. The vertical dashed 
line indicates the position of the PML interface. 

3D Strong blast wave 

Our last test consists of the explosion of a spherical high gas pressure region in a magnetized, initially 
static 3D medium. It has been commonly used for code validation (see, for example Balsa ra fc Spicerl I1999I : 
lLondrillo fc Del Zannal [2000) , and the set-up consists of a cubic domain with 256 grid points in the three spatial 
dimensions spanning from to 1. The initial density, po, is set to unity in all the domain, while the initial pressure is 
set to unity all over except a spheric hot gas region located at the center of the domain of radius ro = 0.125, which is 
a hundred times overpressured (pi = 100). A constant magnetic field with a strength of Box = 10y7^o is initialized 
along the x— direction. 

In Fig. [21] we show a cut of the density in the plane x — y at z = 0.5 and t = 0.02. The system shows the 
axial symmetry imposed by the magnetic field. We can identify the different wave modes present in the simulations. 
The outermost wave corresponds to the fast magnetoacoustic mode, and inside this region there are two wave fronts 
propagating along the magnetic field, which is a slow magnetoacoustic shock. This test verifies that our code can 
handle with three dimensional propagation of highly nonlinear waves at the correct propagation speed, resolving the 
shocks thanks to the hyperdiffusive terms. 
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Fig. 21. — Numerical solution of the 3D strong blast wave density at the plane x — y and z = 0.5 at t = 0.02. From left to right and from 
top to bottom: Mass density, pressure, velocity squared and magnetic pressure. 



